Below you can find a link to the MATLAB code calculating the advection-diffusion problem by the finite element method with stabilization by minimizing the residual for the model Eriksson-Johnson problem.

First, let us recall the advection-diffusion problem, the Eriksson-Johnson model problem (described, for example, in [1] defined on a square domain \( \Omega = [0,1]^2 \) in the following way: We seek \( u \) such as
\( a(u,v)=l(v) \forall v \) where
\( a(u,v) =\int_{\Omega} \beta_x(x,y) \frac{\partial u(x,y) }{\partial x } dxdy + \int_{\Omega} \beta_y(x,y) \frac{\partial u(x,y)}{\partial y } dxdy \\ +\epsilon \int_{\Omega} \frac{\partial u(x,y) }{\partial x} \frac{\partial v(x,y)}{\partial x } dxdy +\epsilon \int_{\Omega} \frac{\partial u(x,y)}{\partial y } \frac{\partial v(x,y)}{\partial y } dxdy \)
\( l(v) = \int_{\partial \Omega } f(x,y) v dxdy \)
\( f(x,y)=sin(\pi y)(1-x) \) is an extension of the Dirichlet boundary condition to the entire area, while
\( \beta = (1,0) \)
represents the wind blowing from left to right, while \( \epsilon = 10^{-2} \) denotes diffusion coefficient. The quantity \( Pe=1/ \epsilon = 100 \) is called the Peclet number, and it defines the sensitivity of the advection-diffusion problem.

In line 830 we specify the diffusion coefficient ϵ=10−2. Its inverse is the Peclet number \( Pe=\frac{1}{\epsilon} \).
We construct the knot vector \( knot\_trial\_x = [0 \quad 0 \quad 0 \quad 1 \quad 2 \quad 2 \quad 2] \) for the approximation space, and vector of point along x axis, between od 0 and 1, \( points\_x = [0 \quad 0.5 \quad 1] \). At these points, we construct the basis of B-spline functions for approximation of the Eriksson-Johnson problem. In a similar way, we define the knot vector for the test space \( knot\_test\_x = [0 \quad 0 \quad 0 \quad 0 \quad 1 \quad 2 \quad 2 \quad 2 \quad 2] \).
According to the idea of the residual minimization method, there must be more testing functions than approximating functions. In our example, we use the quadratic B-spline function with C1 continuity for the approximation space (trial) and the third degree B-spline function with C1 continuity for the test space (test).
The node vectors for the approximating and testing space, and the points on which we span the node vectors, are defined separately for the x axis and for the y axis.

The code execution is also possible in the free Octave environment.
Download code or see Appendix 5A.

The code is activated by opening it in Octave and typing a command
\( advection\_igrm\_adapt \)
After a moment of calculation, the code opens an additional window and draws a numerical and exact solution in it. In the second window, the code draws a residuum, computed additionally by the residuum minimization method. It is a measure of the local solution error.
The code calculates the residuum error \( Residual norm: L2 0.015808, H1 0.100417 \)
The code computes also the L2 and H1 norms of the solution
\( Error: L2 45.90 percent, H1 61.85 percent \)
and compares to L2 and H1 norm of the exact solution.
\( Best possible: L2 2.02 percent, H1 14.20 percent \)
You can see that in order to improve the accuracy of the solution, it is necessary to increase the mesh.